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Classical novae are cataclysmic binary star systems in which the 
matter of a companion star is accreted on a white dwarf (WD) [1, 2]. 
Accumulation of hydrogen in a layer eventually causes a thermonu- 
clear explosion on the surface of the WD [3], brightening the WD to 
~ 10? solar luminosities and triggering ejection of the accumulated 
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matter. They provide extreme conditions required to accelerate par- 
ticles, electrons or protons, to high energies. Here we present the 
detection of gamma rays by the MAGIC telescopes from the 2021 
outburst of RS Ophiuchi (RS Oph), a recurrent nova with a red giant 
(RG) companion, that allowed us, for the first time, to accurately 
characterize the emission from a nova in the 60 GeV to 250 GeV 
energy range. The theoretical interpretation of the combined Fermi- 
LAT and MAGIC data suggests that protons are accelerated to 
hundreds of GeV in the nova shock. Such protons should create bub- 
bles of enhanced Cosmic Ray density, on the order of 10 pc, from 
the recurrent novae. 

A symbiotic nova can be formed when the companion star of the WD is a 
RG. [4]. The ejecta of symbiotic novae expand within the dense wind of the RG 
companion. Novae outbursts usually last from weeks to months. While they 
are expected to repeat hundreds of times [5], the interval between subsequent 
events can be even hundreds of thousand years [6]. However, a subclass of 
objects called Recurrent Novae (RNe) allows one to observe such repeated 
outbursts over a human lifespan [7]. In our Galaxy, ten such objects are known 
in which the repetition of bursts has been seen within a century [6]. According 
to [8] for the symbiotic nova to become recurrent, its WD must be massive 
(21.1 Mo). 

Novae have been deeply studied in the optical and X-ray ranges for decades 
[6, 9-13], but only recently they have been shown as emitters of high-energy 
gamma-ray radiation: first in the case of symbiotic novae [14] and soon after 
with classical novae [15]. Though this clearly indicates that charged particles 
are accelerated to high energies in novae, their nature and radiation mechanism 
are not yet clear. In order to understand the acceleration mechanism of high- 
energy particles, it is crucial to measure the maximum energies of the emitted 
radiation. Until recently, all spectra of gamma-ray novae have been measured 
only up to 6 — 10 GeV range [15] with no hint of emission at higher energies 
[16, 17]. 

RS Oph is a recurrent symbiotic nova with average time between major 
outbursts of 14.7 years [6]. The latest outburst, in August 2021, was promptly 
reported in optical [18] and high-energy (HE, 100 MeV < E < 10 GeV) gamma 
rays by Fermi-LAT [19]. The optical emission showed similar behaviour to 
the 2006 outburst (see Extended Data Figure EDF 1.Following these alerts, 
MAGIC began observations of RS Oph as part of its nova follow-up program 
[17], on August 09, 2021 at 22:27 UT, i.e., about 1 day after the first optical 
and GeV detections. In parallel, the H.E.S.S. collaboration announced very- 
high-energy (VHE, > 100 GeV) gamma rays from RS Oph [20]. The MAGIC 
observations reveal VHE emission contemporaneous to the Fermi-LAT and 
optical maxima, and a decrease below the VHE detection limit two weeks later 
(see Fig. 1). Details of the analysis can be found in Methods section A.1. The 
first four days of MAGIC observations (August 09-12) yield a VHE signal with 
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a significance of 13.20 (see EDF 2), spanning from 60 GeV to 250 GeV, well 
fitted by a single power-law (x?/Naof = 5.9/5). 

Daily spectra are reconstructed (see EDF 3, Method sections A.1 and 
Supplementary section H) allowing us to track the evolution of the outburst. 

The contemporaneous gamma-ray spectrum measured by Fermi-LAT and 
MAGIC can be described as a single, smooth component spanning from 
50 MeV to 250 GeV. Intriguingly, while the GeV emission subsides with a 
halving time scale of — 2.2 days (see also Methods section A.2), the flux mea- 
sured by MAGIC over the first four days is consistent with being constant 
(x2/Naor = 2.9/3), see also EDF 4This suggests a migration of the gamma- 
ray emission towards higher energies, in line with an increase of the maximum 
energies of the parent particles. RS Oph is the gamma-ray nova with the high- 
est flux and energy output to date, as shown by the comparison with the other 
Fermi-LAT detected novae presented in Supplementary section I. Therefore, 
the non-detection of previous novae at VHE range [16, 17] might be explained 
by the lack of sensitivity to dimmer eruptions, without the need to invoke any 
fundamental difference in the spectral energy distribution of RS Oph. 

'The conditions in novae are favourable for the acceleration and subsequent 
emission of radiation by both electrons and protons [15]. The expanding ejecta 
of a nova interacting with the interstellar medium (filled also with the dense 
RG wind in the case of symbiotic novae) will result in the formation of a shock 
wave. Moreover, the fast wind, induced by the nuclear burning on the surface 
of the WD, will catch up with the ejecta, causing an additional internal shock 
[21]. Recently, a correlation between optical and gamma-ray emission has fur- 
ther suggested that a substantial part of the novae explosion's power goes into 
shocks [22]. In such shocks, energetic electrons and protons can be produced 
(see Fig. 2). Gamma-ray emission can arise from photosphere thermal radiation 
up-scattered to the gamma-ray energy range by relativistic electrons via inverse 
Compton scattering. Alternatively, the ambient matter (nova ejecta and RG 
wind) can act as a target for hadronic interaction of protons or Bremsstrahlung 
radiation of electrons [15]. The maximum energies of high-energy particles will 
depend on the efficiency of the acceleration mechanism, duration of the nova, 
and the cooling energy losses (see Methods section B.1 and EDF 5). Protons 
experience only mild cooling by proton-proton interactions with time scale of 
tpp = 21(n,/6 x 105 cm 3) [day], where n, is the number density of the tar- 
get material. Electrons in nova shocks suffer stronger inverse Compton energy 
losses with tro = 4.4 x 10 ?(E/300 GeV)- [1 + 10(E/300 GeV)].? [day]. 
'Therefore, the production of high energy photons via leptonic mechanisms is 
much more demanding on the acceleration processes efficiency than for pro- 
ton models. The simultaneous acceleration of both types of particles (but 
reaching different energies) has also been proposed [17, 23]. We estimate that 
Bremsstrahlung is negligible with respect to inverse Compton component for 
the parameters of RS Oph (see Methods section B). 

We derive the photosphere parameters using fits to the photometry mea- 
surements (see EDF 6) and shock expansion velocity from spectroscopy (see 
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EDF 7). Based on the optical observations of RS Oph during the 2021 out- 
burst, and the derived parameters from previous outbursts of the source, we 
model the gamma-ray emission with the injection of a population of relativis- 
tic electrons or protons (see Methods section B). We take into accout also 
the minor absorption of the emission in the photosphere radiation field (see 
EDF 8). The Fermi-LAT and MAGIC measurement can be well described 
(x?/Naor = 13.1/12, p-value = 0.36) with the proton-only model (see left panel 
of Fig. 3). 

The fit yields a canonical power-law spectrum with an index — —2 and an 
exponential cut-off, corresponding to the maximum energies achieved in the 
acceleration. The day-by-day modeling shows evidence that the energy cut- 
off of protons increases with time (see Supplementary section H and EDF 9). 
'This goes in line with absence of spectral signatures from cooling terms. The 
associated neutrino emission is not expected to be detected by the current 
experiments (see Supplementary section F). 

In contrast, it is difficult to explain the shape of the curvature of the mea- 
sured spectrum between 50 MeV and 250 GeV with leptonic processes. The 
leptonic model requires injection of particles that already contain a strong 
break (change of particles index by 3.25 + 0.28) in the electron energy dis- 
tribution (see Fig. 3, right panel). Since the break must already be present 
in the injection spectrum of particles, it cannot be explained by the cooling. 
In addition, despite a more complicated particle injection model, the descrip- 
tion of the gamma-ray emission in the electron scenario is significantly worse 
(x?/Naor = 27.5/11, p-value = 3.9 x 107° ) than in the case of protons, as can 
be seen in Fig. 3. The relative likelihood of the electron model with respect to 
the proton model for AAIC= 15.3, as defined within the Akaike information 
criterion framework [24], which is normally used for comparison of non-nested 
models, is 4.7 x 1074. 

Despite their intense emission of gamma rays, accelerated protons will even- 
tually escape the nova shock carrying away most of their obtained energy. Such 
protons can contribute to the Galactic Cosmic Rays (CR), which are expected 
to be produced mainly in supernova remnants [25]. 

'The measurement of the proton spectrum required to explain the gamma- 
ray emission of RS Oph can be used to put estimates on novae contribution to 
CR. Using the CR energetics derived for RS Oph (^ 4.4x 109? erg, see Methods 
section B.2), a rate of 50 novae per year [26] would lead to about 0.1% of 
the CR energy contribution from supernovae, which are more rare than novae 
(~ 2 per century) but much more energetic (~ 109? erg). Despite the small 
contribution to the overall CR sea, a nova would significantly increase the CR 
density in its close environment. The energy density of the nova dominates 
over that of the average CR energy density in the Milky Way (^1.8 eV /cm?) 
in a region of radius ~0.5 pc, of the order of the distance to the nearest star 
in our Galaxy. In the special case of recurrent novae, protons accelerated over 
10? yr [27], assuming a recurrent rate of every 15 years, will accumulate in a 
~ 9 pc bubble with enhanced CR density (see Methods section B.2.1). 
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'The detection of gamma rays reaching 250 GeV from a recurrent symbiotic 
nova allowed us to obtain a deep physical insight on the population of rela- 
tivistic particles accelerated by such objects. The modeling of the gamma-ray 
spectrum strongly favors the explanation of the emission via the acceleration 
of protons in a nova shock. Evidence towards the proton acceleration is based 
on: (i) the inferred shape of the energy distribution of injected particles, (ii) 
the better statistical description of the gamma-ray spectral energy distribution 
by the proton model, (iii) the obtained evidence of the increase of the particle 
maximum energies over time, consistent with lack of strong cooling. The pro- 
tons in the nova shock undergo slow cooling, therefore they will be eventually 
able to escape the shock, carrying away a significant fraction of energy. Such 
protons will add to the Galactic cosmic ray budget, however primarily in the 
close neighborhood of novae. 

The observation of the August 2021 outburst of RS Oph introduces a 
new class of sources as VHE gamma-ray emitter: (recurrent symbiotic) novae. 
RS Oph is a recurrent symbiotic nova, the same class of objects as V407 Cyg, 
the first nova detected in the GeV range by Fermi-LAT. While we now know 
that classical novae are also GeV emitters, it is still to be seen if the detec- 
tion of RS Oph emitting in VHE gamma-ray range is due to its recurrent 
symbiotic nature, or just the first sign of such emission from a broader class 
of classical novae. The comparison of gamma-ray measurements in GeV and 
VHE gamma-ray range with previous Fermi-LAT novae does not reveal any 
peculiarity in the emission of RS Oph, except for its brightness (see Fig. 4 and 
EDF 10). Therefore, it is likely that future, more sensitive VHE gamma-ray 
facilities will be able to provide an ample harvest of novae. 


Supplementary information. 
Supplementary sections C-I 
Extended data figures EDF 1-10 
Supplementary Tables 1 — 10 
References (70-97) 
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Methods 


A. Observations and data analysis 


In this section we report the detailed results of the analysis of gamma-ray data 
with MAGIC and Fermi-LAT, and optical data with TJO and ANS. 


A.1 MAGIC 


MAGIC [28] is a stereoscopic system of two imaging atmospheric Cherenkov 
telescopes situated in the Canary island of La Palma, Spain (28.8°N, 17.9°W 
at 2225 m above sea level). Each telescope consists of a 17-m diameter mirror 
dish and a fast imaging camera. The system achieves a sensitivity of (0.92 + 
0.04)96 of the Crab Nebula flux above 210 GeV in 50h in zenith angle range 
30 — 45? [29]. 

MAGIC observed RS Oph in the period between August 09, 2021 to 
September 01, 2021 (MJD 59435.94 to 59458.97) for 34.0 h (see Supplementary 
Table 2). The data quality selection was based on the atmospheric transmis- 
sion and rates of background events. For this analysis we also did not include 
data taken under moonlight condition, as they provide much higher energy 
threshold values. After quality cuts, 21.4 h of the data were used for the analy- 
sis, half of which were taken during the first four days after the nova eruption. 
The source was observed at zenith angles between 36? and 60°. The data were 
taken in the so-called wobble mode, pointing at four different sky positions 
situated 0.4? away from the source to evaluate the background simultaneously. 

The data were analyzed using the MAGIC Analysis and Reconstruction 
Software, MARS [30]. A dedicated low-energy procedure with a special signal 
extraction and image cleaning, the so-called MaTaJu method, was applied 
(see [31] and references therein). Further processing of the data, including the 
image parameterization, the direction and energy reconstruction and gamma- 
hadron separation, were applied following the standard MARS analysis chain. 
The energy threshold of the analysis is —60 GeV. 

We fitted the spectrum obtained from the first four days of observations 
using a single power-law (dN /dE = fo (E/ Eo) °), resulting with a y?/Naor = 
5.9/5 goodness of fit. The used fit also takes into account estimated energy 
bins without detected signal, hence the number of degrees of freedom is larger 
than expected from the number of points in the reconstructed spectrum. The 
normalization energy of the fit (Eo = 130 GeV) is the decorrelation energy (i.e. 
normalization energy which minimizes the correlation of the fit parameters) of 
the four-day sample. The fit parameters are listed in Supplementary Table 3. 

In order to estimate the lower limit on the maximal true energy of gamma 
rays consistent with the MAGIC data we follow the procedure of [32]. We 
perform a likelihood fit of the data with a power-law model with a sharp cut- 
off at a given energy Ec. The 30 (99.7% C.L.) lower limit on the E; is the 
value for which the increase of the x? of the fit is equal to 9. We obtain 170 GeV, 
however taking into account also the 1596 systematic uncertainty on the energy 
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scale following [29] we obtain a slightly less constraining, conservative limit of 
Ecut > 150 GeV. 

We have also performed night-by-night spectral fits to investigate spectral 
variability. The parameters from the first two nights are consistent within 
errors (note however that the exposure on the first night is lower than on the 
remaining ones). A hint of hardening of the emission is seen between the second 
and third night. No significant change of parameters can be seen between the 
third and the fourth night. 

The daily-binned light curve was calculated for an integral flux above 100 
GeV. For the first four days the fit to a constant flux gives a xX?/Nqor = 2.9/3 
with a value of Fç = (4.41 + 0.4644.) x 107! cm? s^ 1. 


A.2 Fermi-LAT 


The Large Area Telescope on-board the Fermi Gamma-ray Space Telescope 
(Fermi-LAT), is a pair conversion telescope designed to detect gamma rays 
with an energy range of 0.02 GeV to > 300 GeV [33]. The Fermi-LAT, with its 
large field of view (2.4sr), observes the entire sky approximately every 3 hours. 
Each analysis is performed with fermitools v2.0.8 and Fermipy v1.0.2 [34] using 
a binned likelihood analysis, P8R3_ V3 instrument response functions (IRFs), 
and the catalog 4FGL-DR2 [35, 36] with the standard Galactic and isotropic 
diffuse background to construct the model of the region of interest (ROI). For 
each analysis, the SOURCE event class is used as this is the recommended 
event class for long duration observations, observations of more than a few 
hours. The SOURCE event class can be further divided into separate event 
types such as PSF0, PSF1, PSF2, and PSF3, where PSF0 corresponds to 
events with the worst PSF and PSF3 are events with the best PSF. 

For the 1-day and 3-day time bins, the Fermi-LAT data-set used encom- 
passes a total time range from MJD 59431.45 to 59461.45, an energy range 
from 0.1 GeV to 1000 GeV, and a 15? ROI centered on the radio coordinates 
of RS Oph (R.A. = 267.555°, Dec. = -6.7078°). We use event type 3, which 
corresponds to all events, for this analysis and select a maximum zenith angle 
of > 90? to reduce any gamma-ray contamination from the Earth limb. The 
majority of 4FGL-DR2 sources for the one and three day time bins are not 
significantly detected (Test Statistic (TS) > 25, see [37]), apart from 4FGL 
J1813.4-1246 and 4FGL J1745.4-0753. These sources correspond to PSR J1813- 
1246 and TXS 1742-078, which are 8.3° and 1.7° away from RS Oph. Here, 
TS is defined as TS = -2 In (Lmaz,0/Lmaz,1), Where Lmaz,o is the maximum 
likelihood of the null hypothesis and Lmaz,ı is the maximum likelihood with 
the source included [37]. The square-root of the TS is approximately equal 
the significance of detection, i.e. a TS = 25 is ~ 50. TXS 1742-078 is a non- 
variable hard blazar and therefore could cause possible source confusion. Due 
to the proximity of 4FGL J1745.4-0753 to RS Oph and possible source con- 
fusion at the lowest energies, the value of the index of 4FGL J1745.4-0753 is 
locked to that of the 4FGL-DR2 catalog. RS Oph is included in the ROI and 
modeled with a Log Parabola model. Additional spectral models were tested 
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for a four-day period contemporaneous to MAGIC observations: a power-law 
(TS= 2168.1) as well as a power-law with an exponential cutoff (TS= 2016.4), 
and the Log Parabola model (TS= 2226.44) had the highest TS, and therefore 
we use the Log Parabola model as our spectral form for RS Oph. The ROI is 
optimized with the normalization and spectral parameters of any 4FGL-DR2 
source with a number of predicted counts < 1 locked to the 4FGL-DR2 values, 
excluding the Galactic and isotropic diffuse background. All parameters on all 
unlocked 4FGL-DR2 sources within 4° are left free to vary, and the ROI is fit 
using Minuit minimizer. If RS Oph source model does not have a TS > 9, num- 
ber of predicted counts > 4 or the error of the integrated flux from 0.1 GeV to 
1000 GeV is greater than 6096 of the value, then it is not considered detected 
and 9596 upper limits (ULs) are calculated. These 1-day and 3-day light curves 
are presented in Supplementary Tables 4 and 5 and in Fig. 1. The 1-day light 
curve in MJD 59435.45-59444.45 range can be well fit (x?/Naot = 6.5/7) with 
an exponential decay with halving time of (2.20 + 0.18) days. 

The analysis of the combined first four days has a data-set which encom- 
passes a time range MJD 59435.45 — 59439.45 and an energy range from 
0.05 GeV to 1000 GeV. Reaching down to 0.05 GeV is necessary to help dis- 
tinguish leptonic and hadronic models described in the main text and seen in 
Fig. 3. The same procedure is applied as in the 1-day and 3-day time bins, with 
some adjustments in the settings to allow the analysis to reach 0.05 GeV. Due 
to the worsening of the Fermi-LAT PSF below 0.1 GeV, we apply a 20? ROI 
centered on RS Oph, and a more restrictive zenith angle selection of > 80?. 
We perform a joint-likelihood analysis with two components, one in the energy 
range between 0.05 GeV to 0.1 GeV and one in the energy range from 0.1 GeV 
to 1000 GeV. We remove PSF0 and PSF1 event types from the analysis below 
0.1 GeV and keep all event types above 0.1 GeV. PSF0 and PSF1 are events 
classified with poor PSF and removing these event types thereby improves the 
PSF with the trade-off of less data. This reduces the possibility of source con- 
fusion from nearby weak sources. This also reduces the chance of false positive 
detections as described in the Fermi-LAT low energy catalog (1FLE)[38]. 


A.3 Optical photometry 


Optical photometric observations of RS Oph were carried out by Joan Oró 
Telescope (TJO) and Asiago Novae & Symbiotic stars Collaboration (ANS, 
telescopes ID 310, 610 and 2203). The TJO is a 1-meter class robotic tele- 
scope located at Montsec observatory (42.05°N, 0.73? E), Catalonia, Spain. The 
multi band (BVR.I.) data were analysed using a semi-automatic pipeline for 
differential photometry [39] assuming the aperture radius of 7.5". The com- 
parison stars magnitudes are obtained from American Association of Variable 
Star Observers International Database (AAVSO). The stars are numbered as 
115, 121, 129, 130, and 133 in the database finding chart. 

The data obtained by ANS are analyzed using PSF photometry method 
described in [40, 41]. The same local photometric sequence, extracted from 
APASS DRS all-sky survey [42, 43] and accurately placed on the system of 
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equatorial standards [44] via the color equations calibrated in [45, 46], has been 
used for all telescopes ensuing a high consistency of the data. The photome- 
try results are given in Supplementary Table6, where the quoted uncertainties 
are the total error, which quadratically combine the measurement error on the 
variable with the error associated to the transformation from the instantaneous 
local photometric system to the standard one (as defined by the photomet- 
ric comparison sequence). All measurements were carried out with aperture 
photometry. 

The cross calibration between instruments was performed by using the 
color index of the source. The data obtained by two telescopes are in good 
agreement. However, to reduce the systematic uncertainties, minimal offsets 
(B —V = +0.03, V — Re = +0.05, and V — I, = —0.02) were applied to 
TJO data. The contribution of the strongest emission lines (H, and Hg) were 
removed from the observed magnitude using the simultaneous spectroscopic 
observations from the publicly available optical spectra in Astronomical Ring 
for Access to Spectroscopy (ARAS) [47]. We found that the contribution of 
the Ha emission line in the V-band is negligible for the first ten days after the 
outburst. The contribution of the H5 emission line is significant in the B-band 
and increases from 396 to 1596 during the same time interval. Moreover, the 
contribution of the H, emission line is dominant in the R-band and increases 
from 596 to 8396 during the same time interval owing to a sudden jump from 
5% to 3496 between T — Ty = 0.98 days and T — Ty = 2.89 days. The results of 
these corrections are presented in Supplementary Table 7. 

All optical data described in this section are corrected for the effect of 
Galactic extinction by assuming E(B — V) = 0.65 [48], Galactic extinction law 
[49], and the absolute fluxes (corresponding to zero magnitude) [50] in each 
band. 

During the nova outburst the photosphere emission creates the dominant 
radiation field. We describe the radiation field using photometric and spec- 
troscopic measurements by applying black body approximation. During the 
first four days of the nova, contemporaneous with the MAGIC measurements, 
the emission can be described by the photosphere temperature dropping from 
Ty, = 10800K to 7680K and radius Ry; = 200 Re (see EDF 6). It should 
be noted that the asymmetry of the photosphere (see e.g. [51, 52], lack of 
measurements at the shortest wavelengths and the presence of lines affect the 
above mentioned fits. Therefore, the photosphere radius and temperature val- 
ues should be considered only a crude approximation of the radiation field, in 
context of gamma-ray emission, and no conclusion on the evolution of those 
two parameters should be drawn. Noteworthly, the photosphere fit of 2006 
eruption [52], when rescaled to the nova distance of 2.45 kpc, provides a similar 
radius (245 — 310) Ro, and temperature (8200 K). 
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A.4 Spectroscopy and ejecta kinematics 


RS Oph spectra during the 2021 outburst have been acquired with the Echelle 
spectrograph of the Varese 0.84 m telescope [53] and the Catania Astrophys- 
ical Observatory Spectropolarimeter [54] of the Catania 0.91 m telescope. 
The reduction of spectra, which included the subtraction of the bias frame, 
trimming, correcting for the flat-field and the scattered light, extraction for 
the orders, and wavelength calibration, was done as in [55] by using the 
NOAO/IRAF packages. IRAF is distributed by the National Optical Astron- 
omy Observatory, which is operated by the Association of Universities for 
Research in Astronomy, Inc. 

The Ha profile obtained on day T — To = 0.91 consists of a triangular shape 
with Full Width at Zero Intensity of ~ 7500 km s^! and a blue shift absorption 
component at 4250kms~!, exactly as it was reported by [56] 1.38 days after 
the 2006 outburst of RS Oph. 

The close similarity of the 2006 and 2021 spectral line profiles along the 
envelope expansion is testified on day T — T + 0 = 15 by the presence of 
satellite components at the same high-velocity (2500 km s^ !). This feature was 
associated by [56] to a presence of two jets (c.f. figs. 1 and 2 therein and EDF 7). 
Also, [57] measured a velocity of 4200 kms! the day after the outburst. 

Because of the day-by-day changing of absorption and emission features 
across the whole RS Oph spectrum, we have determined the velocity of the 
expanding envelope as the terminal value simultaneously representative of the 
Ha, Ha and HeI15876A P-Cygni profiles (Supplementary Table 8). An error 
of 250 kms-! was assumed as representative of differences between profiles. 
EDF 7 shows these profiles in the first three days after the expansion as well 
as on days 5 and 15. 

The acceleration along the initial three days is not statistically confirmed 
and we assume (4500 + 250) kms“! as representative of the ejecta expansion 
at the earliest stage (during the VHE gamma-ray detection by MAGIC). 

It is worth to remind that this velocity is volume average, weighted by the 
brightness, temperature and density of the ejecta velocities and agrees with 
results from the modeling by [58] of the HST images of the spatially resolved 
and expanding ejecta during the 2006 event. Radio maps of the 2006 outburst 
of RS Oph [59] have shown the presence of highly collimated flows with a 
velocity close to 10000 km s^!. In this framework, the decrement of the velocity 
after the initial days is simply a consequence of a non-spherical mass outflow 
[56, 58, 60]. 


B Modeling 


There are compelling both simulation (see e.g. [61]) and observational (see e.g. 
[60]) evidence that the mass transfer in symbiotic binaries causes non-spherical 
circumstellar environment. Such asymmetries are crucial when considering the 
morphology of the emission in particular in optical and X-ray ranges. Here, 
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using a similar approach to [17, 23], we consider a simplified, spherically- 
symmetric scenario in order to evaluate the conditions in which gamma-ray 
radiation can be produced by either electrons or protons and to investigate 
spectral features of such an emission. The used parameters are summarized in 
Supplementary Table 10 


B.1 Acceleration and cooling of particles 


We parametrize the acceleration of charged particles with acceleration param- 


eter €: 
dE || €cE 
(T) Ry " 


where Rr (E) is the Larmor radius of particle with energy E in perpendicular 
magnetic field B. The corresponding acceleration time scale, expressed in days, 
can be computed as: 


tace = E/ (F) =3.9 EJ (Gea) [day]. (2) 


The maximum achieved energies will stem from balancing such acceleration 
time with ballistic time tba, defined as the time from the onset of the nova, or 
by dominating cooling process. The shock distance Rsp at the time t = T — To 
can be estimated based on its speed vgn: 


Ry = 1.2 x 104 ( Sah ) ( t ) [em]. (3) 


4500 km s-1 3d 


As the nova shock expands the adiabatic energy losses will be directly con- 
nected with the ballistic time. We define the adiabatic time scale as the 
time in which the energy of particles decreases by a factor of e, resulting in 
ladiab = €tbal. 

The protons will cool on hadronic interactions with the ambient matter, 
either the nova ejecta, or the RG wind. We assume that the ejecta concentrate 
at the distance of Rsp in a layer with a thickness of h x R,,, with h = 0.1. 
The number density of the ejecta can be estimated as: 


Me; Me; Ush -3 / tN O /hA 1. _ 
J = RL = 60x18 (| 3), 
Pej 4xhR3 m, «10-8 Ma \4500kms—? 3d o1) ol 
(4) 
where Mej is the total ejected mass and m, is the proton mass. Alternative 
assumption that the ejecta fill homogenously a sphere with radius Rsp would 
result in a factor of 3 lower value of nej. The number density of the ambient 
material in the RG wind can be estimated as: 


O Mna 
m 2 
4n Ri, URGMp 


(5) 
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Mna Ush -2 t A URG =1 E 
= 1.1x108 ( ) ( ) 3 
aid 5 x 1077 Mc /yr N4500 km s^! (zz) 1lO0kms-! [em 7], 


where upg is the speed of the RG wind and Mpg is the mass loss rate of the 
RG. The total density of the ambient medium for the hadronic interaction 
for the assumed parameters of RS Oph it is mostly dominated by the ejecta 
(n © Nej). The proton cooling time scale on hadronic p-p interactions can be 
then computed as: 


tpp = (nycagy) ^ = 21(np/6 x 10° cm?) [day], (6) 


where opp = 3 x 1026 cm?. As the cooling timescale is longer than the bal- 
listic time, the maximum energies to which protons can be accelerated are 
determined by the time from the nova onset. 

In the case of electrons, cooling losses can originate either from 
inverse Compton scattering on the photosphere thermal radiation or from 
Bremsstrahlung radiation on the ambient matter. We compute the inverse 
Compton cooling time scale taking into account Klein-Nishina correction factor 
following [62] 


22 2424 1.5 
tro = < (14 den E/ (mec?) ^, (7) 


where me is the electron mass. The total energy density, u>, and characteristic 
temperature of soft photons, ej, can be estimated as 


(Ry, /200 Re)? (Tn /8460 K)4 
(u, /A500 km s—!)2(¢/3 d)? 
2.2(T,n /8460 K) [eV]. (9) 


[erg cm °| (8) 


Uph = 0. 


Eph 


For the used above scaling values the dependence of tro with energy can be 
described as tro = 4.4 x 10 ?(E/300 GeV) !)[1 + 10(£/300 GeV)]-? [day] 
resulting in fast cooling of high-energy electrons. We estimate the 
Bremsstrahlung losses using the same density of ambient matter n, as 


torems = Xo/(NpMpc) = 24(n5/6 x 108 cm~3)~ [day], (10) 


where Xo = 63gcm ? is the radiation length in proton gas. For the expected 
parameters of RS Oph, the Bremsstrahlung losses are thus negligible. Also the 
synchrotron energy losses are negligible, unless the magnetic field in the shock 
reaches the level of about 1G. 

In order to accelerate protons up to energies of a few hundred GeV, the 
value of £B = 1077 G is required (EDF 5). If electrons are accelerated in 
the same conditions, they can reach energies of only ^ 10 GeV. In order to 
explain the observed gamma-ray emission reaching hundreds of GeV, much 
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higher values £B > 3 x 1079 G are required. Second-order Fermi acceleration 
on the nova shock is expected to provide acceleration parameter of the order 
of £ € (vyi/c)? ~ 107+, resulting in the requirement of B > 0.03G fields for 
the electron case and much weaker B > mG for the proton one. 


B.2 Energetics 


The kinetic energy of the ejecta can be estimated as: 


E, — 0.5Me;v2, = 20 x 104 ( Mei ( a ) erg (11) 
m 10-9 M J N4500kms-! 


For the assumed parameters determining the density of target material, the 
fit of the proton energy distribution in Fig. 3 requires a total power in protons 
of 4.4 x 10* erg. This energetics requirement scales with the assumed model 
parameters as: 


aul c Ush 3 d N hA 
Enmon = 044 x 10 (=) Vus me) (zac) Ol 
(12) 
Therefore the efficiency of conversion of energy from the shock to protons can 
be computed as: 


E, nova Me ] =? Ush d =2 h 
= Pro — 09.22 ( — 9 — ( ) 1 
E. Ga (ts ) 4500 km s~! (===) oi 09 


It is clear that protons need to obtain a significant fraction (~ 20%) of the 
shock kinetic energy. Lower fraction could be achieved if the mass of the ejecta 
is higher, it is more concentrated at the shock (lower h) or if the speed of the 
shock is decreased. Concentration of the nova ejecta and proton acceleration 
in the bipolar direction would increase the target material density and effi- 
ciency of the gamma ray production. This would further lower the total energy 
required in the accelerated protons compared to the assumed here spherically 
symmetric scenario. 


B.2.1 Contribution to the Cosmic Ray sea 


These accelerated protons eventually escape the nova to be part of the sea 
of Cosmic Rays. Since they do not suffer strong energy losses due to their 
interaction with intergalactic magnetic and photon fields, as it is the case 
for electrons, their contribution may extend to large distance from the nova 
explosion at all energies. Assuming that the energy released in all novae into 
accelerated protons is similar to that released in RS Oph (Ej, osa = 4.4 x 10? 
erg) and a nova rate of — 50 per year [26] we get a total of 


Novae energy rate = Ep nova x nova rate = 2.2 x 10*°[erg/year] (14) 
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It is considered that a supernova explosion usually releases Esy ~ 10?! erg 
[63], out of which —10% can be converted into accelerated protons at the shock 
between the supernova ejecta and the interstellar medium (ISM). The SN rate 
in the galaxy is ~2 per century [64], therefore the supernova energy rate would 
be: 


Supernovae energy rate = 0.1 x Esn x supernova rate = 2 x 10 5[erg/year] 
(15) 
making the contribution of novae <0.2% to that of supernovae. 

Let us now assume that the average energy density in CRs in the Milky Way 
is Egens,cRs 71.8 eV/cm? [65]. We would like to compute what is the region 
in which the energy density of the protons accelerated by the nova dominates 
over this energy density. The energy density of these protons will be given by 
the total energy (Ep, nova) divided by the volume of the region 


E un 3E, nova 
dens,nova,leruption — 4 Rš 
T eruption 


(16) 


where Reruption is the radius of the region. If we compare Éagensnova = 
Eaens,CRs, We obtain Reruption 0.5 pe, that is subject to the assumption on 
the energy density performed and may change if larger energy densities are 
considered [66]. 

Finally, in the special case of a recurrent nova like RS Oph that repeats its 
explosions every —15 years [67], we would get this energy injection repeated 
over time. Considering a period of recurrence of up to 10? years, the region 
over which this nova would dominate has a size of: 


3Ep,nova X 104 


ËEjdens,nova,recurrent = 4 R3 
T Itrecurrent 


(17) 


and the radius over which the protons accelerated by the nova would dominate 
over the energy density of the ISM would be Ryecurrent ~9 pc. 


Availability of data and materials: Analysis products of MAGIC data are 
available here: http:/ /vobs.magic.pic.es/fits/. Low level data are available on 
request. 

Code availability: The code for fitting the electron and proton models is 
available in https:/ /opendata.magic.pic.es/download?pid—2. 
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Fig. 1 Multiwavelength light curve in the VHE (MAGIC, top panel), high-energy (Fermi- 
LAT, middle panel) and optical (TJO, ANS, and AAVSO, bottom panel) bands. The lack 
of MAGIC data between MJD 59440 and MJD 59454 is due to the presence of bad weather 
conditions and strong moonlight. Errorbars represent 1-sigma statistical uncertainties in the 


data points. 
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Fig. 2 Schematic representation of RS Oph during an outburst. A photosphere (yellow 
circle) surrounds the White Dwarf (WD, white small circle). Its companion star, a red giant 
(RG, red circle) emits a slow wind (red arrows). Ejecta of the nova explosion (gray arrows) 
propagate into the surrounding medium causing a shock wave encompassing the binary 
system (gray dashed line). In the shock wave, energetic electrons and protons (magenta and 
green wavy lines, respectively) are trapped by a magnetic field and accelerated. Gamma 
rays (white arrows) are produced by either electrons scattering the thermal radiation of the 
photosphere (yellow arrow) or by protons interacting with the surrounding matter (gray and 
red dots). 
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Fig. 3 Gamma-ray spectrum of RS Oph observed with Fermi-LAT (empty crosses) and 
MAGIC (filled circles) averaged over the first four days of the outburst, modeled within 
hadronic (left panel) or leptonic (right panel) scenario. The dashed line shows the gamma 
rays from the 7° decay and the dotted line shows the inverse Compton contribution of 
the secondary e^ pairs produced in hadronic interactions. dN/d E, and dN/d Ee report the 
shape of the proton and electron energy distributions obtained from the fit. The bottom 
panel shows the fit residuals. Errorbars represent 1-sigma statistical uncertainties in the data 
points. 
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Fig. 4 Total energy vs duration of RS Oph 2021 outburst compared to that of the other 
novae detected by Fermi-LAT. Data taken from [14, 15, 68]. Errorbars represent 1-sigma 
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Extended Data Figure EDF 1 Optical B-band observed magnitude and the color index 
of RS Oph 2021 outburst from ANS (red circles) and TJO (green triangle) compared to that 
of 2006 eruption (black crosses, computed with respect to MJD of 53775.86, [69]). Top three 
panels show the color indices, while the bottom panel shows B magnitude evolution. 
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Time - 10.80 h 

Nan = 3975; N,, = 3085.4 + 31.7 
Nex = 889.6 + 70.6 

Significance (Li&Ma) = 13.20 
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Extended Data Figure EDF 2 Distribution of the squared angular distance between the 
nominal source position and the reconstructed arrival direction of events (black crosses) and 
the estimated background (gray shaded area). Vertical dashed line represents the angular 
cut below which the number of background and excess events as well as the statistical 
significance of the detection are given (inset panel). Error bars represent 1-sigma statistical 
uncertainties in the data points. 
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Extended Data Figure EDF 3 Modeling of daily emission in proton model (left panels) 
and electron model (right panels) for first, second, third and fourth day after the nova 
eruption (from top to bottom). The dashed line shows the gamma rays from the 7° decay and 
the dotted line shows the inverse Compton contribution of the secondary e* pairs produced 
in hadronic interactions. dN/d E, and dN/dE¢ report the shape of the proton and electron 
energy distributions obtained from the fit. The bottom panel shows the fit residuals. Error 
bars represent l-sigma statistical uncertainties in the data points. 
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Extended Data Figure EDF 4 Comparison of the photon flux measured by Fermi-LAT 
above 100 MeV with the one measured by MAGIC above 100 GeV (left panel) and with that 
of the V-band obtained by ANS (right panel). Arrows show the sequence of the flux temporal 
evolution and the blue line shows the linear proportionality fit in both panels. Error bars 
represent 1-sigma statistical uncertainties in the data points. 


Springer Nature 2021 IATEX template 


Proton acceleration in thermonuclear nova explosions revealed by gamma rays 37 
protons, £B=1.0e-07 G electrons, £B-1.0e-07 G 
19] ballistic 10? 
=== acceleration 
- adiabatic : 
- pp 10: 
3 a 
z = sl Pur vm ass t q PET 
$ a 
cul Aces oT Mese ec) SMO! EEE Y Sus oru mui e 
2 2i 
B^ d eneaudsvayduissai wna Gas tiade ues s 
5 8 
3 E 
» 10? 29 10? 
: Ei. 4 o = ballistic 
s acceleration 
10-4 107 =-=- adiabatic 
—-- Bremsstrahlung 
— IG 
107? 107 
10? 101 102 103 104 105 10? 10? 10? 10? 10^ 105 
Energy [GeV] Energy [GeV] 
protons, £B=3.0e-06 G electrons, €B=3.0e-06 G 
101 mE mer ——————— s eee 
Z 10° É 10° 
$ a 
3 3 
4 E: 
3 107 G 107 
3 T 
H £ 
š E o o fe ballistic 
= 40-2 F ig 
10737 7 0 y 0 [e ballistic 10 s acceleration 
== acceleration --- adiabatic 
dos --- adiabatic ia —-- Bremsstrahlung 
==. pp — c 
10? 10? 10? 103 104 105 10? 101 10? 10? 10^ 105 
Energy [GeV] Energy [GeV] 


Extended Data Figure EDF 5 Cooling and acceleration time scale for protons (left 
panels) and electrons (right panels) for two values of £ B parameter: 1077 G (top panels) and 
3 x 10-9 G (bottom panels). Assumed parameters (see text for details): vs, = 4500 kms-!, 
t = 3d, Rpn = 200 Ro, Tjj = 8460 K, np = 6 x 108 cm-?. 
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Extended Data Figure EDF 6 Optical photometry performed by ANS 1, 3, and 4 days 
(see the panel titles) after the outburst (blue empty markers) corrected for the Galactic 
absorption. Filled markers show the flux after subtraction of Ha and Hf line contributions. 
'The thick black lines show a black-body emission used in the modeling, while the dashed line 
shows for comparison the average 2006 spectral fit from [52] (with the photosphere radius 
corrected to the distance of 2.45 kpc). Horizontal error bars represent the bandwidth of the 
filters used. 
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Extended Data Figure EDF 7 Example of Ha, Hg and HeI5876A P-Cygni profiles used 
to determine the behavior of the expansion velocity of the expanding envelope (time after 
the outburst is given in the top right part of each panel). The bottom right panel shows the 
evolution of the velocity in time. Error bars represent 1-sigma statistical uncertainties in the 
data points. 
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Extended Data Figure EDF 8 Absorption of the gamma-ray emission on the radiation 
field of the photosphere and collision with it. Assumed parameters: vs, = 4500kms-!, 
Ryn = 200 Re. Temperature of the photosphere is T,, = 10780 K, 9490 K, 8460 K and 
7680 K for the time after the nova onset: 1d (black solid), 2d (red dotted), 3d (green dashed), 
4d (blue dot-dashed) respectively. 
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Extended Data Figure EDF 9 The maximum energy of protons obtained from the 


theoretical model fits to the daily gamma-ray emission 


points) shown in Fig. 3. Red and 


green line show, respectively, the scenario of proportional increase and constant value of 


maximum energy. Error bars represent 1-sigma statistical 
of the maximum energy of protons. 


uncertainties in the determination 
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Extended Data Figure EDF 10 Comparison of RS Oph to other Fermi-LAT-detected 
novae. Spectra of other Fermi-LAT-detected novae are shown in the top panel. Gamma-ray 
spectra of V407 Cyg (middle panel) and V339 Del (bottom panel) compared to the measured 
(red) and scaled (gray) RS Oph gamma-ray spectra. Blue triangles and arrows correspond 
to Fermi-LAT measurements and upper limits of V407 Cyg (top) and V339 Del (bottom). 
Red squares are the Fermi-LAT spectrum of RS Oph and red circles the MAGIC one. Gray 
squares are the Fermi-LAT scaled spectrum of RS Oph and gray circles the MAGIC one. 
Cyan arrows correspond to the VERITAS (V407 Cyg) and MAGIC (V339 Del) upper limits. 
'The dashed blue lines correspond to the best-fit using a single power-law for the Fermi-LAT 
data. 'The dotted blue lines correspond to the best-fit using a power-law with an exponential 
cut-off for the Fermi-LAT data. Data taken from [14-17, 68]. Error bars represent 1-sigma 
statistical uncertainties in the data points. 


Springer Nature 2021 IATEX template 


Proton acceleration in thermonuclear nova explosions revealed by gamma rays 
Supplementary material 


C The recurrent symbiotic nova RS Oph 


RS Oph is composed by a massive carbon-oxygen white dwarf (WD) [70] and 
a M0-2 III mass-donor RG star [71]. The orbital solution implies a WD mass 
of Mw p = 1.2 — 1.4Mc and an RG mass of Mra = 0.68 — 0.80Mo. This nova 
has shown eight eruptions between 1898 and 2006 [6]. Interestingly, RS Oph 
was pointed out as a plausible source from which GeV emission can be detected 
[72]. [73] showed that the prompt hard X/soft gamma-ray emission of the 2006 
outburst of RS Oph detected by Swift/BAT [11] could not be accounted by the 
decay of radioactive isotopes. [72] proved that this emission could be explained 
via the production of non-thermal particles by diffuse shock acceleration and 
proposed RS Oph as a possible GeV candidate. It is also a type Ia Supernova 
progenitor candidate [6, 74]. The system has a period of (453.6 0.4) days [75]. 
The system has a circular (e z 0) orbit[76, 77], however a mild eccentricity 
(e = 0.14 + 0.03) has been claimed as well [75] probably due to a better cover- 
age of the radial velocity curve. The estimation of the wind mass loss rate of 
the RG in RS Oph is ~ 5 x 1077 Mc yr-! [61]. The total matter ejected dur- 
ing nova outburst is difficult to estimate. Models of the explosion give values 
2 x 1077 — 1078 Mọ for different values of the WD mass [78]. A rough esti- 
mation of “not much more than 1077 Mọ” has been given by [10], however it 
was based on earlier measurements of the RG wind density that is an order of 
magnitude lower than the one in [61]. 

An outburst was reported on August 08, 2021 (MJD To — 59434.93) [19]. 
'The spectroscopy measurements performed in the first few days of the nova 
show mild acceleration of the ejecta from 3700 — 2700 km s^! at Ty + 0.87 d to 
4200 — 4700 km s^! at To +2 d for the Ha and H8 P Cyg lines, respectively [79, 
80]. We assume an ejecta speed of vsn = 4500 km s^! , see Methods section A.4. 
The equipartition magnetic field derived from observations starting 20 days 
after the 2006 outburst is B = 0.08 — 0.11 G [81], while the estimate 18.5d 
after the 1985 nova onset was 0.04 G [82]. In the case of a similar recurrent 
nova, V1535 Sco, a value of B = 0.13 — 0.17 G was measured one week after 
the outburst [83] and for V745 Sco 0.03 G at the distance of the shock of 
4.5 x 1014 cm [84]. It should be noted that, as the shock dissipates, B declines 
over time. 


C.1 Estimates of the RS Oph distance to Earth 


The distance to RS Oph has been object of intense debate (see Supple- 
mentary Table 1). Historically, a value of 1.6 kpc was estimated [85, 86] and 
canonically assumed. Dedicated discussions on the distance were performed in 
the past [87] in which a distance of 1.4 kpc was considered the most likely one. 
The value is however at odds with the mass accumulation rate needed for rep- 
etition period of RS Oph. Namely, at such assumed distance, the calculated 
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blackbody radius of the secondary star must greatly underfill its Roche lobe 
[88]. 

More recently, the parallax distance to the source of (2.68 + 0.16) kpc was 
provided by Gaia [89]. However, as it is discussed in [88], Gaia Data Release 2 
(and Gaia Early Data Release 3) do not have reliable measures of the parallaxes 
for RS Oph. The issue arises due to the long-period binary orbit which makes 
the center of light wobble back and forth with a greater amplitude than the 
parallax itself, not providing a good fit to the single-star model applied in 
these data releases. 

Using the argument that the RG needs to fill its Roche lobe to efficiently 
accrete matter onto the WD, the favored distance to the source is 3.1 + 0.5 kpc 
[87]. [67] also pointed out that using the light curve information, the most 
likely distance is 4.3 + 0.7 kpc, and the earlier, lower estimates suffered from 
overestimated absorption along the line of sight (see the discussion in [88]). 
'This is however at odds with the expansion velocity of the synchrotron shock, 
as it was pointed out by [90], in which they derived a distance of (2.45 + 
0.37) kpc. Given all these caveats, in the subsequent calculations we assume 
the distance to be 2.45 kpc. 


D Multiwavelength view 


In EDF 4, we compare the integral fluxes (obtained from daily spectral fits) 
of Fermi-LAT and MAGIC. Fitting the relation with a linear proportionality 
(F(» 100GeV) < F(> 0.1GeV)) we obtain x?/Naor = 12.7/3. Therefore 
simple, achromatic gamma-ray variability is unlikely (chance probability p — 
5.3 x 1073). 

We also perform a joint fit of Fermi-LAT and MAGIC data with a 
Log Parabola function shape: dN/dE = fo x (E/Eg) % ñln(É/Eo) with 
Eo —130 GeV. The fits are summarized in Supplementary Table 9. 

In the right panel of EDF 4, we compare the integral fluxes (obtained 
from daily spectral fits) of Fermi-LAT and the differential fluxes for the V- 
band obtained by ANS and TJO. We selected only the coincident data for 
MJD 59435 - 59445, for which the daily Fermi-LAT integral fluxes fulfill the 
condition not to calculate an upper limit. We fit a function y = Az to the 
data that gives a xX?/Nqor = 9.8/7. 

Despite the data being consistent with a linear optical-GeV correlation, 
this does not provide a straight-forward interpretation about the underlying 
particle population producing this emission. As the ejecta propagate away 
from the WD, the radiation field seen at the shock will in fact decay faster 
than the observed at Earth optical emission. Moreover, as the electrons would 
cool faster (see Methods section B) than the observed decay of the radiation, 
we do not expect a linear relation between the IC emission and optical. The 
IC emission would rather follow the rate of injection of particles from the 
acceleration process. This correlation could therefore be related to ejecta or 
particle acceleration parameters. 


Springer Nature 2021 IATEX template 


Proton acceleration in thermonuclear nova explosions revealed by gamma rays 


— IceCube limits 
-8 -vV 


e 


Vu 


-8.5 


TTTITTTTTTTTTTTTTTTTTTTT[TTTTRITTTTTTTTTTTTT 


-1 0 1 2 3 4 


5 6 7 
log (Energy / GeV) 


Supplementary Figure 1 Predicted neutrino emission (ve with cyan dashed line and v,, 
with solid blue line) associated to the proton model (see Fig. 3) compared with 9096 C.L. 
limits obtained by the IceCube Collaboration [91]. 


E Absorption of gamma-ray radiation 


In both electron and proton models, the production of the gamma-ray radiation 
occurs relatively close to the photosphere Rs = 2.8(t/[day]) x Rpn, a strong 
thermal source. We investigate the effect of absorption of the produced gamma 
rays on such a radiation field. We derive angle-dependent optical depths in 
vicinity of such thermal source and compute the average absorption as 


Absorption(E) — os f e T (E89) sin 5, d6,j,, (18) 
0 


where 05; is the angle between direction of photon and the radial direction from 
the center of the photosphere. Additionally we assume that photons crossing 
the photosphere are fully absorbed (7 = oc). 

Derived absorption of the emission is presented in EDF 8 for different 
days after the nova outburst. In general the absorption is not very strong, in 
particular at energies S 300 GeV, where gamma-ray emission was detected. 
Nevertheless it is taken into account in the modeling. 


F Neutrino emission 


The gamma-ray emission in hadronic scenario would be accompanied by neu- 
trinos. We calculated the neutrino emission corresponding to the proton model 
presented in Fig. 3 and compared it with limits from the IceCube Collabo- 
ration [91]. It is clear that due to sub-TeV energies achieved by protons, the 
predicted neutrino emission does not reach energies higher than those of pro- 
tons and these limits cannot constrain the model (see Supplementary Fig. 1. 
We also investigated if SuperKamiokande could have detected neutrino emis- 
sion associated to the nova outburst. However, due to low collection area at 
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Supplementary Figure 2 Fit to the Fermi-LAT and MAGIC SED with a proton-electron 
model. The dashed line shows the gamma rays from the x? decay and the dotted line shows 
the inverse Compton contribution of the secondary e= pairs produced in hadronic interac- 
tions. dN/dE, and dN/dE¢ report the shape of the proton and electron energy distributions 
obtained from the fit. The bottom panel shows the fit residuals. Errorbars represent 1-sigma 
statistical uncertainties in the data points. 


the GeV energies [92] the expected number of events is only of the order of 
5 x 1077. 


G Proton-lepton model 


The presence of high-energy protons or electrons is not only dependent on 
their maximum energies (see Methods section B.1). Differences in the injection 
process of electrons and hadrons into the acceleration mechanism (see [93] 
and references therein) can cause preferential dominant acceleration of one 
or the other type of particles. Following |17, 23], we test as well a model 
in which both electrons and protons are accelerated in the same shock. We 
assume injection with a power-law and exponential cut-off for both particles 
types. The cut-off energies are related by the cooling/acceleration balance (see 
Methods section B.1). The resulting best fit is presented in Supplementary 
Fig. 2. The assumed spectral shape of injected electron and proton populations 
cannot explain the emission well (x?/Naor = 29.4/11, corresponding to p-value 
of 2.0 x 10 3). The best fit also requires L, / L, ~ 2, much larger than < 0.1 
constrained in observations of V337 Del [17]. 
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H Day-by-day proton modeling 


In addition to the modeling of an average state of the source in the first 4 
days, we also perform modeling of individual days after the nova onset. The 
results of fits with the proton model are shown in left panels of Fig. 3. On 
individual days the preference of the proton model over the electron model is 
lower, however except for the first night, the electron model provides lower p 
value. Summing up over the four days x2 increases by 5.2 despite additional 4 
parameters, results in AAIC= 13.2 which corresponds to AIC likelihood ratio 
of 1.4 x 1073. 

Interestingly, the spectra show a hint of gradual softening of the power-law 
component accompanied by and increase of the value of the cut-off energy (see 
Fig. 9). Such behaviour is in line with the expectations from the cooling and 
acceleration time scales defined in Methods section B.1. Namely, due to low 
cooling losses of protons, their maximum energies are mainly determined by 
the duration of the acceleration. The dependence of the maximum energy of 
protons on time can be fit very well (x2 /Nao£ = 0.54/3) with such a scenario of 
proportional increase with time (corresponding to £B = const). Such continu- 
ous increase of maximum proton energies could last until the shock is drained 
up from its energy, or is slowed down by the interstellar medium. However, 
as the target material dilutes with time, the expected gamma-ray emission 
would fall below the detectability level. Constant value of the cut energy can 
be excluded at chance probability 3.1 x 10^? level (x?/Naor = 13.9/3). It 
should be noted that while the fit only considers statistical uncertainties of 
the reconstructed maximum energy, it is unlikely that any systematic uncer- 
tainties would mimick such a gradient as the data are taken over a time span 
of only a few days in similar observational conditions. 


I RS Oph in context with other novae 


I1 Gamma-ray novae 


To put the RS Oph eruption into context, we compared it to other published 
Fermi-LAT detections of novae: V407 Cyg 2010 [14], V3124 Sco 2012, V959 
Mon 2012, V339 Del 2013 [15], V1369 Cen 2013 and V5668 Sgr 2015 [68]. There 
are other studies of Fermi-LAT novae [22], apart from several ATels and sub- 
threshold sources [94] of classical and symbiotic novae that are not included 
in the comparison presented in this section. It is important to mention that 
although RS Oph is considered to be a symbiotic nova, it was pointed out [95] 
that even though RNe have nova eruptions on symbiotic stars, they may not 
share all the properties of symbiotic novae (in particular very slow and low 
amplitude eruptions without Roche lobe overflow). 

On the top panel of EDF 10, we present a comparison of the RS Oph Fermi- 
LAT SED coincident with the MAGIC detection (MJD 59435.8 - 59439.8) 
and the average of the full flare (MJD 59434.8 - 59464.8) compared to the 
aforementioned novae. We can see that both the flux corresponding to the 
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simultaneous data, and the average flux during the whole eruption are from a 
factor of a few up to almost two orders of magnitude larger than previously- 
detected eruptions. 

To perform meaningful comparisons, we defined the duration of RS Oph 
eruption determined by the intervals spanned by the TS > 4 in the daily light 
curves [68]. With this definition, the duration is 30 days, that is comparable 
to the rest of the Fermi-LAT published novae. This has not only been the 
eruption with the highest flux, but also the most luminous one as it is shown 
in Fig. 4, for which we have used the results of the fit with an Exponential 
Cut-Off Power Law fit for E > 0.1 GeV of the average flux from the eruption 
as in [68]. This statement is dependent on the assumed distance of 2.45 kpc, 
and is subject to the uncertainties in the determination of the distances to 
different novae (see the discussion in Methods section C.1). 

In [68] there is the speculation that there is an apparent inverse relationship 
of the total energy with gamma-ray durations for classical novae, that would 
also roughly be followed by V407 Cyg. The fact that we measured that RS Oph 
has a factor of a few higher energy emitted in gamma rays, points to intrinsic 
differences between this eruption and the others detected in classical or symbi- 
otic novae. The total power of gamma rays emitted from RS Oph, 1.8x10?? erg 
that is about 0.9 x 107?(M.;/(10~° Msun)7 * (Usp /4500kms 1) ?(d/2.45 kpc)? 
of the kinetic energy of the shock. 


I.2 Detectability of novae at VHE gamma-ray range 


We perform a comparison of the spectrum of RS Oph eruption with the most 
similar nova detected at gamma rays so far: V407 Cyg [14]. In the top panel 
of EDF 10, we can see the comparison between the average V407 Cyg spec- 
trum measured by Fermi-LAT during the 22 days of its eruption and the ULs 
by VERITAS [16] for a total of —5 hour observation time, compared to the 
MAGIC and the Fermi-LAT flux simultaneous to the MAGIC detection of 
RS Oph. We also scaled MAGIC and Fermi-LAT RS Oph flux to reach that 
of V407 Cyg measured by Fermi-LAT. We can see that in every case, the 
UL established by VERITAS on V407 Cyg lies above the extrapolation of the 
RS Oph flux measured by MAGIC, therefore the RS Oph results are in agree- 
ment with the non-detection by VERITAS, assuming that the VHE gamma-ray 
emission from V407 Cyg follows the same spectral shape as that of RS Oph. 
There are physical differences between classical novae and RS Oph, a 
recurrent nova with a strong wind from the companion, that could cause the 
difference in gamma-ray emission. To evaluate the detectability of classical 
novae, we nevertheless performed a comparison between RS Oph and V339 Del 
[17], observed by MAGIC during its eruption. In the bottom panel of EDF 10, 
we can see the average V339 Del spectrum measured by Fermi-LAT during 
the 27 days of the eruption and the MAGIC ULs compared to the RS Oph 
measurement of MAGIC and Fermi-LAT simultaneous to that of MAGIC. We 
also scaled the RS Oph flux to reach that of V339 Del measured by Fermi- 
LAT for the simultaneous data to MAGIC. We can see that the MAGIC ULs 
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for V339 Del are below the RS Oph measurement, however, if we scale RS Oph 
flux down, the MAGIC ULs of V339 Del are above the MAGIC measurement 
of RS Oph. We note the caveat of comparing the average fluxes measured by 
Fermi-LAT for previous novae and that simultaneous to MAGIC measurement 
for RS Oph, in which the state was high. 

We can conclude that the detection of RS Oph at VHE gamma rays was 
possible due to a higher gamma-ray flux, rather than favorable spectral dis- 
tribution shape, and that other previously detected novae could have emitted 
photons up to the same energies, that remained undetected due to the sen- 
sitivity of the observations. This means that a VHE gamma-ray instrument 
more sensitive in the — 100 GeV energies would open the possibility of detect- 
ing a large number of gamma-ray emitting novae if their emission extends up 
to VHE [25, 96]. 
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Supplementary Table 1 Different distances estimated for RS Oph. The distance 
assumed in subsequent calculations is marked with asterisk. 


Distance [kpc] Method Reference 
1.6 H I absorption measurements [85, 86] 
luos Several estimations [87] 
2.45 + 0.37* Expansion velocity [90] 
3.1 + 0.5 Requirement of RG filling its Roche lobe [87] 
4.3 + 0.7 Light curve [67] 
2.68 + 0.16 Parallax [97] 


Supplementary Table 2 Summary of the MAGIC observation campaign: time of 
observation slot, observation conditions, total observation time during the slot, effective 
time after the data selection (only Dark data). 


MJD Start - End 


Time after cuts |h] 


59435.94 - 59435.98 
59436.89 - 59437.04 
59437.89 - 59438.03 
59438.88 - 59439.02 
59439.89 - 59440.02 
59440.89 - 59441.02 
59444.89 - 59444.91 
59445.88 - 59445.90 
59451.89 - 59452.00 
59452.88 - 59453.01 
59453.88 - 59454.00 
59454.87 - 59454.98 
59455.87 - 59455.97 
59456.87 - 59456.97 
59458.89 - 59458.97 


Obs. conditions Obs. time |h] 
Dark 1.0 
Dark 3.6 
Dark 3.2 
Dark 3.2 

Dark + Moon 3.0 

Dark + Moon 3.0 
Moon 0.1 
Moon 0.2 

Dark + Moon 2.1 

Dark + Moon 2.9 

Dark + Moon 2 
Dark 2.5 
Dark 2.8 
Dark 2.3 
Dark 1.9 


1.0 
3.5 
3.1 
3.2 


Supplementary Table 3 Spectral fit results (normalization fo at Eg = 130 GeV, and 
photon index a) of MAGIC data: daily and the combined emission from the first four days. 
Errors represent 1-sigma statistical uncertainties in the fits. 


| MJD fo [10-19 TeV-1cm-7s-1] a X^ /Naot 
59435.94 - 59435.98 5.011 3.927078  5.6/5 
59436.89 - 59437.04 3.135 oh 4.717035 51/5 
59437.89 - 59438.03 5.0377 84 3.70025 -316/5 
59438.88 - 59439.02 4.83017 3.781025 — 10.3/5 
59435.94 - 59439.02 456 0 4.07103 — 59/5 
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Supplementary Table 4 Fermi-LAT 1-Day Average Integral Flux. Values with * are 
fixed in the UL calculation. Errors represent 1-sigma statistical uncertainties in the fits. 


Integral Flux (> 0.1 GeV) 


MJD Start - End TS (10-7 photeusem 2s-1) a B 
59431.45 - 59432.45 0.4 « 11.0 2.0* 0.0* 
59432.45 - 59433.45 0.0 « 10.8 2.0* 0.0* 
59433.45 - 59434.45 0.0 « 10.7 2.0* 0.0* 
59434.45 - 59435.45 191.2 18.8 + 3.1 2.16 + 0.13 0.054 + 0.076 
59435.45 - 59436.45 — 1006.9 46.4 + 3.9 1.96 + 0.078 0.197 + 0.051 
59436.45 - 59437.45 501.0 39.3 + 4.3 2.123 + 0.099 0.175 + 0.066 
59437.45 - 59438.45 — 433.4 27.4 t 3.4 1.955 + 0.095 0.169 + 0.065 
59438.45 - 59439.45 197.8 17.9 + 3.4 2.12 + 0.16 0.24 + 0.12 
59439.45 - 59440.45 172.3 12.8 + 2.7 1.96 + 0.16 0.22 + 0.11 
59440.45 - 59441.45 94.6 7.8 + 2.3 1.63 + 0.21 0.16 + 0.11 
59441.45 - 59442.45 97.6 12.2 + 3.2 1.99 + 0.17 0.15 + 0.12 
59442.45 - 59443.45 60.3 4.3 + 1.6 1.96 + 0.42 0.95 + 0.58 
59443.45 - 59444.45 82.0 5.5 + 1.9 1.58 + 0.28 0.39 + 0.21 
59444.45 - 59445.45 2.5 « 11.4 2.0* 0.0* 
59445.45 - 59446.45 24.9 10.0 + 4.1 2.65 + 0.61 0.18 + 0.33 
59446.45 - 59447.45 16.0 « 11.9 2.0* 0.0* 
59447.45 - 59448.45 17.1 < 11.7 2.0* 0.0* 
59448.45 - 59449.45 28.5 < 12.1 2.0* 0.0* 
59449.45 - 59450.45 3.2 < 11.5 2.0* 0.0* 
59450.45 - 59451.45 27.6 4.0 + 2.3 1.93 + 0.34 0.16 + 0.22 
59451.45 - 59452.45 0.7 < 11.2 2.0* 0.0* 
59452.45 - 59453.45 2.4 < 11.1 2.0* 0.0* 
59453.45 - 59454.45 13.3 < 11.3 2.0* 0.0* 
59454.45 - 59455.45 20.2 6.1 + 3.1 2.48 + 0.49 0.15 + 0.28 
59455.45 - 59456.45 5.0 < 11.6 2.0* 0.0* 
59456.45 - 59457.45 7.5 < 11.7 2.0* 0.0* 
59457.45 - 59458.45 1.0 « 11.0 2.0* 0.0* 
59458.45 - 59459.45 5.8 UT 2.0* 0.0* 
59459.45 - 59460.45 6.1 « 11.6 2.0* 0.0* 
59460.45 - 59461.45 12.8 4.2 t 2.3 2.26 + 0.46 0.51 + 0.26 
59461.45 - 59462.45 8.3 « 11.3 2.0* 0.0* 
59462.45 - 59463.45 4.4 « 11.4 2.0* 0.0* 
59463.45 - 59464.45 4.6 « 11.5 2.0* 0.0* 
59464.45 - 59465.45 1.9 « 10.8 2.0* 0.0* 


Supplementary Table 5 Fermi-LAT 3-Day Average Integral Flux. Values with * are 
fixed in the UL calculation. Errors represent 1-sigma statistical uncertainties in the fits. 


Integral Flux (> 0.1 GeV) 


MJD Start - End TS (10-7 photonscm- 28-1] a B 
59431.45 - 59434.45 0.0 « 10.5 2.0* 0.0* 
59434.45 - 59437.45 1621.7 34.0 + 2.1 2.041 + 0.053 0.160+ 0.040 
59437.45 - 59440.45 797.5 19.7 + 1.8 1.999 + 0.072 0.200 0.051 
59440.45 - 59443.45 — 244.0 8.2 + 1.5 1.83 + 0.12 0.200 0.084 
59443.45 - 59446.45 80.6 4.0 + 1.3 1.90 + 0.22 0.28+ 0.16 
59446.45 - 59449.45 49.0 1.62 + 0.67 1.47 + 0.40 0.55 0.32 
59449.45 - 59452.45 25.4 2.3 + 1.3 2.01 + 0.32 0.16 0.20 
59452.45 - 59455.45 27.4 2.0 + 1.2 1.93 + 0.29 0.03 0.13 
59455.45 - 59458.45 11.3 « 10.8 2.0* 0.0* 
59458.45 - 59461.45 10.8 « 11.6 2.0* 0.0* 
59461.45 - 59464.45 9.0 « 11.2 2.0* 0.0* 
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Supplementary Table 6 The observed optical magnitude of RS Oph. Errors represent 
l-sigma statistical uncertainties in the measurements. 

MJD Telescope B V Re Ic 
59435.913 ANS 5.667 + 0.011 4.884 + 0.014 4.194 + 0.016 3.544 + 0.018 
59437.824 ANS 6.611 + 0.011 5.855 + 0.014 4.816 + 0.016 4.164 + 0.020 
59437.850 ANS 6.508 + 0.013 5.811 + 0.022 4.801 + 0.028 4.167 + 0.032 
59438.820 ANS 6.968 + 0.034 6.321 + 0.045 4.966 + 0.046 4.611 + 0.054 
59438.845 ANS 6.871 + 0.012 6.215 + 0.024 4.982 + 0.036 4.459 + 0.030 
59439.812 ANS 7.288 + 0.014 6.637 + 0.017 5.270 + 0.019 4.856 + 0.020 
59439.881 TJO 7.303 + 0.012 6.605 + 0.011 - 4.850 + 0.015 
59440.817 TJO - 6.895 + 0.015 5.482 + 0.014 5.103 + 0.011 
59440.818 ANS 7.555 + 0.008 6.923 + 0.009 5.421 + 0.011 5.050 + 0.014 
59440.840 ANS 7.445 + 0.008 6.808 + 0.014 5.285 + 0.026 4.911 + 0.021 
59441.817 TJO - 7.137 + 0.020 5.587 + 0.014 5.291 + 0.013 
59442.886 TJO - 1.328 + 0.011 - - 
59442.936 TJO 7.974 + 0.012 7.335 + 0.011 5.704 + 0.011 5.431 + 0.011 
59443.824 ANS 8.053 + 0.008 7.436 + 0.010 5.745 + 0.011 5.537 + 0.014 
59443.886 TJO 8.014 + 0.012 7.359 + 0.011 5.740 + 0.011 5.537 + 0.011 
59444.313 ANS 8.129 + 0.008 7.508 + 0.009 5.797 + 0.011 5.618 + 0.014 
59444.926 TJO 8.184 + 0.013 7.616 + 0.011 5.881 + 0.011 5.719 + 0.011 
59445.849 TJO - 7.688 + 0.012 5.967 + 0.011 5.782 + 0.011 
59446.826 ANS 8.367 + 0.010 7.752 + 0.012 6.003 + 0.014 5.875 + 0.016 
59446.897 TJO 8.381 + 0.013 7.695 + 0.011 6.011 + 0.011 5.867 + 0.011 
59447.807 ANS 8.582 + 0.008 7.893 + 0.010 6.115 + 0.012 6.037 + 0.014 
59447.917 TJO 8.557 + 0.014 7.914 + 0.012 6.105 + 0.011 6.024 + 0.011 
59449.808 ANS 8.655 + 0.008 7.961 + 0.014 6.157 + 0.016 6.139 + 0.014 
59449.899 TJO 8.719 + 0.028 8.193 + 0.019 6.255 + 0.013 - 
59452.824 ANS 8.993 + 0.007 8.431 + 0.009 6.374 + 0.014 6.492 + 0.015 
59454.882 TJO 9.027 + 0.012 8.433 + 0.011 - - 
59456.813 ANS 9.106 + 0.008 8.532 + 0.016 6.645 + 0.017 6.728 + 0.018 
59458.804 ANS 9.181 + 0.016 | 8.594 + 0.019 6.535 + 0.048 6.774 + 0.038 
59459.801 ANS 9.277 + 0.009 8.748 + 0.010 6.792 + 0.011 6.984 + 0.012 
59459.805 ANS 9.305 + 0.016 — 8.695 + 0.018 6.641 + 0.051 6.899 + 0.041 
59459.834 ANS 9.385 + 0.007 8.879 + 0.008 6.862 + 0.009 7.043 + 0.014 
59460.801 ANS 9.403 + 0.009 8.848 + 0.015 6.894 + 0.016 7.086 + 0.017 
59460.808 ANS 9.461 + 0.007 8.958 + 0.008 6.913 + 0.010 7.148 + 0.014 
59460.871 TJO 9.361 + 0.013 8.918 + 0.0011 6.972 + 0.011 7.134 + 0.011 
59461.825 ANS 9.597 + 0.007 8.997 + 0.014 7.039 + 0.015 7.230 + 0.016 
59462.796 ANS 9.690 + 0.033 9.070 + 0.043 6.995 + 0.070 7.228 + 0.054 
59463.805 ANS 9.627 + 0.008 9.087 + 0.011 7.095 + 0.012 7.307 + 0.012 
59463.883 ANS 9.789 + 0.010 9.296 + 0.018 7.299 + 0.021 7.459 + 0.024 
59465.791 ANS 9.781 + 0.017 9.161 + 0.021 7.089 + 0.034 7.363 + 0.039 
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Supplementary Table 7 The observed optical magnitude (B and Re band) of RS Oph 


after removing the contribution of Ha and Hg emission lines. Errors represent 1-sigma 
statistical uncertainties in the measurements. 

MJD Telescope | Hg/ B (96) B H. Rec (96) Re 
59435.913 ANS 3 5.698 + 0.011 5 4.247 + 0.016 
59437.824 ANS 9 6.701 + 0.011 31 5.134 + 0.016 
59437.850 ANS 9 6.598 + 0.013 34 5.093 + 0.028 
59438.820 ANS 11 7.081 + 0.034 48 5.393 + 0.046 
59438.845 ANS 11 6.984 + 0.012 48 5.409 + 0.036 
59439.812 ANS 11 7.405 + 0.014 48 5.694 + 0.019 
59439.881 TJO 11 7.420 + 0.012 - - 
59440.817 TJO - - 69 6.054 + 0.014 
59440.818 ANS 12 7.676 + 0.008 69 5.993 + 0.011 
59440.840 ANS 12 7.566 + 0.008 69 5.857 + 0.026 
59441.817 TJO - - 87 6.268 + 0.014 
59442.936 TJO 13 8.103 + 0.012 91 6.407 + 0.011 
59443.824 ANS 14 8.192 + 0.008 89 6.434 + 0.011 
59443.886 TJO 14 8.153 + 0.012 89 6.429 + 0.011 
59444.313 ANS 14 8.268 + 0.008 89 6.486 + 0.011 
59444.926 TJO 14 8.330 + 0.013 83 6.538 + 0.011 


Supplementary Table 8 Log-book of spectroscopic observations, and expansion velocity 
of the Ha, Hg and HeI15876A P-Cygni profiles. A conservative error of 250 km s-! has 
been associated to all velocities. 


Telescope Serra la Nave Varese 
Spectrograph CAOS Echelle 
R-—A/AA 45 000 18 000 
Range 400-900 nm 425-890 nm 
MJD Expansion velocity [km s `1] 
59435.837 4250 
59436.820 4600 
59437.807 4750 
59438.830 4000 
59439.808 3000 
59440.867 2800 
59441.810 2700 
59442.838 59442.824 2700 
59443.821 59443.806 2700 
59444.850 59444.810 2600 
59445.852 59445.853 2500 
59446.817 59446.808 2500 
59447.801 2500 
59448.796 2400 
59449.794 2400 
59450.822 59450.814 2400 
59451.796 2400 
59452.785 2400 
59454.804 2300 
59455.792 2300 
59459.858 2100 
59467.830 2100 
59470.835 2000 
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Supplementary Table 9 Daily Fermi-LAT and MAGIC joint spectral fit results. 


Individual columns give normalization fo at normalization energy Eo = 130 GeV, slope a 
at Eo, curvature parameter and goodness of fit (x? /Naor). Errors represent 1-sigma 


statistical uncertainties in the fits. 


| MJD follo 1 TeV-Tcm-7s-1] a B x7/Naot 
59435.94 - 59435.98 5.4 1.3 3.86 + 0.13 0.194 + 0.019 6.1/6 
59436.89 - 59437.04 4.54 + 0.78 3.73 £0.11 0.175 + 0.020 16.4/6 
59437.89 - 59438.03 5.37 + 0.85 3.64 + 0.12 0.173 + 0.020 3.7/6 
59438.88 - 59439.02 5.00 + 0.78 3.44 + 0.14 0.147 + 0.027 10.8/6 

| 59435.94 - 59439.02 5.08 + 0.45 3.697 + 0.059 0.175 0.010 9.3/6 | 


Supplementary Table 10 Summary of the nova parameters used for the modeling of 
nova gamma-ray four day averaged spectrum (see Fig. 3). Parameters marked with an 


asterisk have modified values in the night-by-night modeling (see EDF 3) 


Parameter Symbol Value 
Distance d 2.45 kpc 
Photosphere radius Rph 200 Ro 
Photosphere temperature Tph 8460 K* 
time after nova explosion t 3d* 
Expansion velocity Ush 4500 km s-! 
Mass of nova ejecta Mej 1076 Mo 
Confinement factor h 0.1 


